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In this paper we propose a perturbative method for the reconstruction of the covariance 
matrix of a multinormal distribution, under the assumption that the only available information 
^ amounts to the covariance matrix of a spherically truncated counterpart of the same distribution. 

I— ' We expand the relevant equations up to the fourth perturbative order and discuss the analytic 

^ properties of the first few perturbative terms. We finally compare the proposed approach with an 

^ exact iterative algorithm (presented in ref. p[I) in the hypothesis that the spherically truncated 

\^ covariance matrix is estimated from samples of various sizes. 

(N 

1 Introduction 

O 

CN In a recent paper [1], we have studied how the covariance matrix {&B)ij = co\{Xi, Xj \ X G B^{p)) 

of a multinormal random vector X = {Vfc}^^-|^ ~ A/'t,(0, S) in i; > 1 dimensions, conditioned to a 
centered spherical domain By{p) = {x £ W" : x^x < p}, relates to the unconditioned covariance 

^ matrix Sjj = cov(Xj, Vj). Thanks to the symmetries of the geometrical set-up, and S can 

^ be shown to commute. Moreover, if we denote respectively by // = {/ifc}fc=i and A = {Afc}^^]^ the 

eigenvalue spectra of and S, then 

fJ'k = M—, — TV' k = l,...,v, (1.1) 

a{p;X) 

with a and belonging to the class of Gaussian integrals 

r 2 2 2 ^ — y^/(2r)) 

aMm..Xp-A)= / d^x^^^...n'^(^.,A,), 6{y,v) = ] . (1.2) 



'Corresponding author, e-mail: filippo.palombi@istat.it 
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Since &s and S are simultaneously diagonalizable, we can assume E = diag(A) with no loss 



of generality. Reconstructing S from ©g means solving eq. (1.1) with respect to A under the 
assumption that p and fi are given. This is only possible provided /i fulfills specific algebraic 
constraints (see sect. 2 of |.lj). Anyway, a closed-form solution is out of reach, due to the non- 
linear character of the problem. For this reason, we have proposed in [1] a numerical technique 
based on a fixed point iteration, whose convergence mechanism is related to some conjectured 
correlation inequalities within Bv{p)- 

In this paper, we approach eq. ( |1.1[ ) from a different perspective. We move from the observation 
that a simplified set-up occurs when the eigenvalue spectra are fully degenerate, a case which has 
been first considered by Tallis j2j. If ^ui = . . . = /x„ = /2, by symmetry it follows Ai = . . . = At, = A. 



The converse holds true as well. Eq. (1.1) reduces in this limit to 



A = A^(f) -W), (1.3) 

with (x) denoting the c. d.f. of a x^^variable with v degrees of freedonj^ It can be readily checked 
that the function Tp{X) is monotonic increasing in A. In addition, we have 

i) lim Tp{X) = , ii) lim TpCX) = ^— , (1.4) 

A^O A-s>oo V + I 



whence we recognize that eq. (1.3) can be numerically inverted (by any root-finding algorithm) 
provided < /i < p/(w + 2). 



No w, e q. (1.3) can be thought of as the lowest order approximation of a perturbative expansion 
of eq. ( 1.1 ) around the point At = {A, . . . , A}. If the condition number of S is not extremely large, 
such an expansion is expected to quickly converge, so that a few perturbative corrections to At 
should be sufficient to guarantee a good level of approximation. 

In critical applications, requiring reconstructions of the covariance matrix for several values of 
p or /i, it could be important to access fast yet approximate solutions, such as perturbation theory 
provides, rather than slow yet exact ones. Indeed, the convergence of the fixed point iteration has 
been shown to slow down as p — )• or f — )• oo, the rate of the slowing down being polynomial in the 
former limit and exponential in the latter. Thus, in all situations where p <^ min^jAfc} oi v ^ 1, 
the use of the fixed point algorithm could be unfavorable. 

Another advantage of the perturbative approach turns up when is not known exactly, 
but instead it comes along as the result of a multivariate Gaussian sampling from a spherically 
truncated population of finite size. In that case, we shall see that the statistical fluctuations of 



the higher components of p are amplified by eq. (1.1) as a consequence of the non-linearity of the 



problem, sometimes resulting in unacceptable variances for the higher components of A. In the 
framework of perturbation theory, non-linearity arises systematically on increasing the order of the 
approximation, since each perturbative correction depends non-linearly upon the previous ones. 
Therefore, by stopping the expansion at different orders, we have the possibility to define a class 
of statistical estimators of A, each characterized by its own bias and variance. The variance of the 
upper (lower) half of the spectrum increases (decreases) along with the perturbative order, whereas 
the bias always decreases. In all applications where the the upper part of the spectrum counts, it 



^with abuse of notation we shall often write %r^(f ) in place of _F„+2(p/A)/_F„(p/A). 
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is therefore possible — at least in principle — to optimize the choice of the perturbative estimator 
according to one's needs. 

Aim of the present paper is to carry out a theoretical study of the perturbative expansion of 



eq. (1.1) up to the fourth order and to discuss some analytic aspects of it. Here is a plan of the 
paper. In sect. 2 we derive some preliminary results concerning the integrals a and a^, which are 
necessary for a systematic implementation of the perturbative strategy to all orders. In sect. 3 
we work out the expansion by means of paper-and-pencil calculations and Maple"^'^ programs 
(given in the appendix). Sect. 4 is devoted to discussing some analytic properties of the first few 
perturbative terms, when /i is chosen to be the average of the truncated eigenvalues. Finally, in 
sect. 5 we simulate in sample space the statistical properties of the perturbative estimators of A and 
the iterative one for a specific choice of the covariance matrix. Conclusions are drawn in sect. 6. 



2 Building blocks in Tallis' limit 

As widely known, perturbation theory is an expansion technique around a reference solution, which 
is assumed to be either calculable or easily computable. Its mathematical structure develops from 
building blocks which are themselves entirely defined in terms of the reference solution. When 



it comes to perturbing eq. (1.1) around At, the elementary objects we need to focus on are the 
Gaussian integrals ak£m... and their partial derivatives in the limit A — t- At, which we shall also refer 
to in the sequel as Tallis' limit. 

To begin with, let us set-up the notation. Throughout the paper we shall keep on represent- 
ing a fully degenerate spectrum by At = {A, . . . , A}. Outside Tallis' limit an appropriate index 
reshuffling allows us to assume with no loss of generality the ordering < H2 ^ ■ ■ ■ ^ IJ-v (which 
is induced by Ai < A2 < . • . < A^,, as asserted by Proposition 2.3 of [1]). We shall denote a 
generic Gaussian integral with n (not necessarily distinct) indices {ii, . . . = I either by the 
notation aij^,,,i„{p', A) introduced in sect. 1 or by the compact notation (p; A), where sub- 

script colons are meant to separate each directional index from the multiplicity it has in I. By 
the same token, we shall denote the n*^ order derivative operator with respect to the variances 
Ail ... , -^in either by the standard symbol = d'^/{dXi^ . . . dXi^) or by its compact version 

9i:mi...«:m„ = a™i+-+'^V(5Ai)'"i . . . (^A^)™-. Of course, the multiplicity set Mx = {mi, . . . ,m4 
in unambiguously associated to I. For consistency reasons it must fulfill X]fc=i "^fe = ^- "^^^e 
of vanishing multiplicities we shall drop all the corresponding indices. For instance, we shall 
write afc:mfc(p;A) in place of the rather pedantic ai:o,,,k;mk...v:o{p'i ^) as well as dk-.nn. in place of 
di:o...k:m^...v:0- Last but uot least, we refer the reader to [1] for definitions and properties concerning 
the truncation operator Tp and its inverse r~^. 

Having said that, we start our investigation of Tallis' limit with a simple proposition, which is 
just meant to review the findings of [2]: 

Proposition 2.1. For mi > 0, . . . , m^, > 0, we have 

'^l:mi...v:mv{p] Xt) = ^l:mi...v:mv Fv+2{mi+...+mv) (t^ ' (2-1) 
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with 



Ai:™,...„:„„ = J|(2mj-1)!!. (2.2) 



Proof. In order to derive eq. (2.1), we simply need to represent ai:mi...i):m„ (/o; '^t) in spherical 
coordinates, 

xi = rfi{9i, . . .,9v-i) , 

; (2.3) 
Xv = rfv^Oi, . . .,6v-i) ■ 

Recall that Ylk fk~^ ^^'^ '^^■^ ~ '^''"^drdO, with dJl embodying the angular part of the Jacobian 



of eq. (2.3) and the differentials of the angles 9i, . . . , A few plain algebraic passages yield 



«l:mi...'[;:m„(p; At) — Ai:mj...t,:m^F„_l_2(mi + ...+m„) ) (2-4) 

with the proportionality coefficient A.i-mi...v:my being independent of p or A. In order to fix it, we 
observe that ai:-mi...t):m„(/o; At) factorizes in v one-dimensional integrals as p — >• oo, corresponding 
to unconditioned univariate Gaussian moments of orders 2mi,. . . ,2mt,, normalized respectively by 
powers A™^,. . . ,A'"" of the common variance. Hence, we infer A.i:mi...v:my = 11^=1(2"^*; — l)-- Q 

The x^-c.d.f. F^+2(mi+...+mi,)(p/A) is intuitively interpreted as a correction factor incorporating 
all the effects of the volume conditioning. Obviously, it lessens the value of the integral except when 
p — 7- oo, where it becomes ineffective. 



Whenever the Gaussian integral is expressed in standard notation, eq. (2.1 ) can be only applied 
provided the index multiplicities are preliminarily counted. Nevertheless, it is reasonable to expect 
Ai:mi...i):m„ to be the compact representation of some not yet specified coefficient Ajj,,.j„. By the 
same argument as above this is recognized to be 

Ai„„i„=E[zl...zl], provided M{0,1) , k = l,...,n. (2.5) 



Thanks to Isserlis' Theorenjj [3j , Aj^...j^ can be reduced to a sum of products of Kronecker symbols, 
namely 

A,,..,„=Y,ll^ij, (2.6) 

where the sum extends over all distinct ways of partitioning the duplicated set 2^ = {ii, ii, . . . ,in, in} 
into pairs. For later convenience it is worthwhile listing a few frequently recurring examples: 

An = 1 , (2.7) 

A,^i, = 1 + 26i,i, , (2.8) 



Eqs. (2.5)-(2.6) allow us to reformulate Proposition 2.1 according to 



^in the framework of mathematical physics the same result is universally known as Wick's Theorem. 
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Corollary 2.1. For n >0 and {ii, . . . a set of n (not necessarily distinct) indices, we have 



a 



A 



ii...i„-P^v+2n 



(2.10) 



with A 



ll...lr, 



as in eq. (2.6). 



Now, in order to take arbitrarily high order derivatives of aj^,..j„ (p; A) with respect to the 
covariance eigenvalues, we can iterate the basic rule 



(2Afc5fc) an...i„(/3; A) = aii...i„fc(p; " + 2 '^fei,^ "n. ..»„(/>; A) , 



(2.11) 



which follows by differentiation under the integral sign. Note that eq. (2.11) is not specifically 
related to spherical truncations, i.e. it is formally invariant under a reshaping of the truncation 
surface. It also shows that the differential operator 2X^8^ behaves in a simpler manner than 
when acting on aij^„,i^^, in that it produces an integer linear combination of similar integrals. For 
this reason, the recurrence generated by d^-^n can be derived by first working out the action of the 
operator {2Xkdk)"' and then using 



dk-.r 



2i 



(2.12) 



where the symbols ["] denote unsigned Stirling numbers of the first kind. Eq. (2.12) is a classic 
of combinatorial analysis. The reader is referred for instance to exercise 13, chap, b of ref. [1] for 
a proof of it. Iterated applications of the operator 2X^3^ generate increasingly involved sums of 
Gaussian integrals, as asserted by 



Proposition 2.2. For all j > 1 and n> 0, we have 



{2Xkdkyak;n = ^(-1)^ ''Cjr{n)ak:{n+r) , 



r=0 



with 



r Ix 



'^i-r-\ J-r 



(2.13) 



(2.14) 



Proof. The proof is by induction. We first note that cio(n.) = (2n + 1) and cii(n) = 1. Hence, 
for j = 1 eq. (2.13) agrees with eq. (2.11). Now, suppose that {2Xkdky Ok-.n is well represented by 
eq. (2.13) with Cjr{n) as in eq. (2.14). Then, 



{2Xkdky~^^ak;n = '^{-'^y ''Cjrin){2Xkdk)ak:(n+r) 
r=0 
j 

= ^{-'^y~''cjr{n){ak:(n+r+i) " [2(n + r) + l]afc.(„+^)} 

r=0 

i+1 

= Y.^-'^y'^^~'{cj(r-i){n) + [2(n + r) + l]cj>(n)}afc,(„+,) 



(2.15) 



r=0 



5 



Hence, the proof is complete if we are able to show that Cj>(re) fulfills the recurrence 

C(j+i)r{n) = Cj(,,_i)(n) + [2(n + r) + l]cj>(n) , 
To this aim, we first calculate the second term on the r.h.s. as 

r £i ^j-r-i j—r 

[2(n + r) + l]c,,(n) = ^ E ' ' • E [^(^ + + 1)] U^^'^in + 4) + 1] 

= EE---E [2(n + r) + l)] n [2(^ + 4) + 1] 



(2.16) 



£2=0^3=0 £j+i_^=0 



s=2 



= EE--- E n [2(-+4) + i], 

£l=r£2=0 £j+i_r=0 s=l 

and then we add it to the first one, thus obtaining 

r-l ei t-j-r j+l-r 

c,-(,_i)(n) + [2{n + r) + l]c,.(n) = E E ' • ' E 11 [^(^ + + 1] 

^1=0^2=0 ej + l-r=0 S = l 

+ E E • • • E n + ^s) + 1] = C(,+i),(n) . 



(2.17) 



(2.18) 
□ 



Nested sums similar to eq. (2.14) are considered for instance in [5j, where all the cases in 
study are reduced to closed-form expressions with the help of special numbers, such as binomial 



coefficients, Stirling numbers, center factorial numbers, etc. Eq. (2.14) looks a bit harder to manage, 



since the summand is a product of non-homogeneous functions of the sum variables, whence it is 
not clear whether the nested sum can be ultimately evaluated in closed-form. However, as far as 
we are concerned, a perfectly convenient representation of Cjv(n) is provided by 

Proposition 2.3. For all j > 1, < r < j and n > 0, we have 



Cjr[n) 



[2(n + t)-l]!! (j 



t=r 



[2(n + r) - 1]!! [t 
with the symbols {]} denoting Stirling numbers of the second kind. 



(2.19) 



Proof. Let us denote by djr{n) the r.h.s. of eq. (2.19). For j = 1, we have dio{n) = 2n + 1 



and dii{n) = 1, which equal respectively cio(n) and cii(n). Thus, we just need to prove that 



djr{n) obeys eq. ( |2.16 ). To this end, it is sufficient to make use of the basic recursive formulae 
{"^'} = ^C} + lii} and ("+^) = a + L"_i). We detail the algebra just for the sake of 
completeness. We start from 



, oy+l-f [2(n + t)-l]!! + (t 



t=r 
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3-t 



t=r 



[2(n + t)-l]!! ^.Aj\(t 
[2{n + r)-l]\V '\ti\r 



+ E (-2) 



[2(n + i) + l]!! fiUt + l 



t=T-l 



[2(n + r) - 1]!! \t 



{2n + l)djr{n)+ (-2) 



[2(n + t) + l] !! fj) ^ t 



t=r-l 



[2{n + r) - 1]!! \t] \r - 1 



(2.20) 



Then, we add and subtract {2r)djr{n) to the r.h.s. of eq. (2.20). Hence, 
<i(j+i)r.(n) = [2{n + r) + l]djr{n) 



j-t 



[2(n + t)-l]!! (j 



t=r-l 



[2{n + r) - 1]!! [t 



[2(n + t) + l] 



r — 1 



2r 



(2.21) 



Since r(*) = (t — r + 1)(^,^-|^), the second term on the r.h.s. is recognized to be (ij(j._i)(n). □ 



In view of [5] it is no surprise that Cjr{n) can be represented in terms of Stir. 



ing numbers. In 



We recah indeed that Stirhng numbers of the first and second kind fulfill the identity 

max{j,fc} 

, I r I K' 



addition, the presence of terms such as {-J} fits perfectly when combining eq. (2.12) with eq. (2.13). 
dl indeed that Sti: 

E } 

t=o ^-^ ■' 



(2.22) 



We have thereby collected all the ingredients needed to prove the main result of this section, 
namely 

Proposition 2.4. For m,n > 0, let fC = {ki, . . . , km} and X = {ii, . . . denote sets of (not 
necessarily distinct) indices, i.e. 1 < kj < v and 1 < ij < v. Then, 



1 



9fei...fc„aii...i„(/3; At) ^ 

In particular, for n = 0, 1 we have 

1 



-A 



J7T,?1 ...271 



j=0 ^ ^ 



dki...kmO^{p; K) 
5fei...fe„ai(p; At) 



2™ A" 
1 

2m 



-A, 



ki...k 

n 



E(-:r-<7).«.(e), 



-A 



g'-'"i7)^— (!) 



(2.23) 



(2.24) 



(2.25) 



Proof. First of all, let Aij^ = {mi, . . . , my] and M.x = {ni, . . . , n„} be the multiplicity sets associ- 
ated respectively to K. and X, such that dk^...k^ = di:mi...v:my, ah...i„ = ai:ni. and 

dki...km'^h...in{P', At) = 5i:mi...-u:m„ai:ni...i);n„(p; At) • (2.26) 
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Recall also that M.]c and M.x fulfill Yl\=i — and X^^^^ iT-e = n- for consistency reasons. Using 
eq. (2.12) m times yields 



5fci...fc„an...i„(/3; At) = r- X] 



ii=i 



2Ji 



mi 



- (-1) 



2J« 



ruv 



■.ni...v:ny 



(2.27) 



Moreover, with the help of eqs. (2.13) and ( |2.l[ ), eq. (2.27) reduces to 
'9fei...fc„aii...i„(/3; At) 



I mi -I r -1 Jl rriv ^ 

J_ J_ "ii f.iymi-gi _L 

2^1 71 ^ • ■ ■ 2> 

ii=l 91=0 jr„=l 

[2(ni + gi)-l]!! ... [2K + Q,)-1]!! 



gt,=0 



(2.28) 



As a next step, we evaluate the coefficients Cj^q^(ni), Cj^q^{ny) in terms of the expressions 
obtained in Proposition 2.3. Then, we make use of eq. (2.22) to make all Stirling numbers disappear. 



We finally identify 11^=1 [2("i + mj) - 1)]!! = /^i:{rni+ni)...v:{m,+n^) and so arrive at 



9fei...fe„an...i„(/3; At) 



2™ A" 



-A 



l:(mi+ni)...f :(m„+n„) 



mi m„ 



91=0 <ji,=0 



j-(gi+...+g„)/"il 
91 



t)+2(s+(ji+...+(ji,) , , 



(2.29) 



Two additional moves are still needed to complete the proof. In first place, we notice that the 
multiplicity set associated to /C U X = {ki, . . . , km, h, . . . , in} is Mkux = {"ii + ni, . . . , rriy + n^,}. 
Therefore, owing to Isserlis' Theorem we can identify '^i:(m+mi)...v:{n^+-m^) = ^k^...k^h...in,- In 
second place, we observe that ((?i + . . .+qv) ranges from to m when < gi < mi, . . . , < gt, < m^. 



Accordingly, we can recast the r.h.s. of eq. (|2.29|) to 
5fci...fc^aii...i„(/3; At) = 



1 / P A 

^Afc,,..fc„ii..,i„ ^(-l)'""^ejF^+2(s+j) (^-j 



(2.30) 



with 



mi m„ 

E-E 

gi=0 g„=0 



mi 



m, 



(2.31) 



This multiple sum can be easily calculated by reiterating Vandermonde's convolution (^) (p^^) = 
("■+") , which finally yields = C^) . □ 
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3 Perturbative expansion 



In order to invert eq. (1.1), perturbation theory prescribes that we interpret fj, and A as smooth 
functions of a parameter e G [0,1], such that A(e = 0) = At and A(e = 1) = • fi. We must 
consider e as an auxihary variable, allowing us to pass continuously from Tallis' limit to the ultimate 
solution we are looking for. Then, we are supposed to expand //(e) and A(e) in power series of e 
around the point e = 0, namely 



A + 6A« + e2A(^)+... 



(3.1) 



/xfc(e) = /i + e//^')+eVf^+... . (3.2) 
For later convenience we shall denote by Rk the integral ratio a^/a. Since depends smoothly 



upon A, it can be analogously expanded in power series of e. Thus, eq. (1.1) reads 



R + eR^K^Rf + ... 



(3.3) 



with R = Rk{p', At) and R 



(n) 



"^d"i?fc/de"|e=o for n = 1, 2, . . . The idea underlying perturba- 



tion theory is that we treat separately terms belonging to different perturbative orders, that is to 



say we equal terms of the same order in e on both sides of eq. (3.3) and then solve one by one the 
algebraic equations obtained. 

There are of course some caveats. 

i) Since is an input parameter, an implementation of the perturbative strategy is only possible 
provided we prescribe how it enters the Taylor coefficients of the function //^(e). In principle, the 
assignment can be decided in complete freedom. For instance, the choice which will be made in the 
sequel is to concentrate on the lowest available Taylor coefficient, i.e. 



I Mr 



0, n = 2,3, , 



(3.4) 



Alternatively, we could spread /j.^ over all Taylor coefficients, e.g. according to the prescription 



(0) 
(n) 



1 

nl 



{iog[i + (/.fc-/i)]r 



n 



1,2,... 



(3.5) 



Both eqs. (3.4) and (3.5) comply with the requirement /Ufc(e = 1) = fik- Nevertheless, it must be 
borne in mind that each legitimate splitting of ^ affects differently the convergence properties of 
the perturbative series of A as well as the statistical properties of the truncated series, when is 
turned into a random variable in sample space. This will be further investigated in sect. 5. 

ii) The specific choice of /i is rather arbitrary: as far as we are concerned with the feasibility of 



the perturbative expansion, the only requirement to fulfill is that eq. (1.3) be invertible, which is 
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guaranteed provided < fi < p/{v + 2). A convenient choice, as we shall see in the sequel, is 
represented by 



1 



(3.6) 



k=l 



Proving that /i < p/{v + 2) whenever /i G ^{Tp^) is not trivial. Such an upper bound relies indeed 
on Gaussian correlation inequalities similar to those conjectured in sect. 2 of [T]. Here we give 
a tentative proof. If ^ G V{t~^), then it follows fik = \k{oik/oi) for some A E W" . In order to 
establish an upper bound to /2, we have to inspect how this varies as a function of A, so we look at 
the derivatives 



1 " 

dkfi = - dkPi 



i=l 



dkPk + ^ dkPi 



dk {^k^) + ^>'idk 



a 



2vX 



^ cov {Xi,Xi I X G B,{p)) , (3.7) 



where the rightmost expression has been obtained with the help of eqs. (2.9)-(2.11) of ref. [Ij. In 
that paper we conjecturecj^ cov | X € Bv{p)) < for A; 7^ i. Here, we invoke the additional 

conjecture 



var {Xi I X G B,{p)) > ^ |cov {Xf,Xl \ X G B,{p)) \ , 



which we motivate by the intuitive observation (supported with no exception by extensive numerical 
tests) that the negative square correlations induced by conditioning X to B.^ (p) are extremely weak. 



From eqs. (3.7)-(3.8) it follows 



y^Uk< lim V'Afc — (/>;A) 

^ (Ai,...,A„)^(oo,...,oo)^ a 



vp 
v + 2 



(3.9) 



This estimate improves the one given in eq. (2.17) i) of ^ and consequently leads us to conjecture 
V{t~^) C H'^{p) C Hy{p), where Hy{p) has been defined in [1] and 



H^{p) = jxGM!; : Xfc <min|^,— ^1 and < ^| 

is a tighter bounding region. 

iii) Perturbation theory works only provided the 0(e") -equations 



(3.10) 



j=0 



n = 0, 1, . . . 



(3.11) 



^this conjecture makes only sense under the assumption X ~ JVv{0, E) with E = diag(A). 
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generated by collecting all the 0(e")-terms from eq. (3.3), establish an algebraic relation among 



the Taylor coefficients of A(e) which can be solved with respect to A^"^. This allows us to repre^ 
sent the latter as a function A('^''(A, A^^-*, . . . , X^'^~^\ n^""^) of its lower order companions and 



If such property is confirmed, as we shall argue in a moment, then inverting eq. (1.1) pertur 



(IT 



fc=i' 



batively amounts to solving one by one in sequence the systems of equations 

(2) 

... up to a predefined order n. Establishing the level of precision thus achieved is 
a complicated matter, as typical when facing perturbative expansions. From a qualitative point of 
view, the approximation is certainly recognized to be correct up to 0(e"'''^)-terms, which however 
is not an estimate of the truncation error. 

With regard to point Hi), we notice that the only contributions to if^"^ holding terms which 



are proportional to A^"^ are precisely A^" i? and . All the other contributions, of the form 

n — 1, can only depend upon A, A(i), A("-i). Indeed, depends 



A 



("-i) 7?(i) 



RX' with j = 1, . . . 



upon e implicitly via A(e), thus its j*^ order derivative with respect to e distributes progressively 
according to the chain rule of differentiation. When evaluating such derivative at e = 0, all terms 
proportional to strictly positive powers of e vanish. As a consequence, every surviving term must 
be proportional to a product of Taylor coefficients of A(e), each belonging to {A, A*^^), . . . , A^-'^}. In 
particular, when j = n an explicit calculation yields 



terms in A^'^) within : R + ~^Y1 >^f^djRk[p; At) . 



(3.12) 



In order to evaluate the first order partial derivatives of Rk, we make use of Propositions 2.1 and 2.4, 



djRkip; Xi 



a 



1 

2A 



(1 + 25 J k) 



i?2 



26 



'jk- 



+2 



F, 



(3.13) 



whence eq. (3.11) reduces tcr 



(n) 



(a,a«,...,a(-i),/.(-)) , 



(3.14) 



with 



•3^ kj 



(1 + 25k j) 



t)+4 



fVh2 

F2 



(3.15) 



We have thereby obtained a system of linear equations with A*-"^ and Q^'^'> representing respec- 
tively the variable and constant vectors. Moreover, the coefficient matrix J is the Jacobian of the 
truncation operator Tp in Tallis' limit. Its determinant is given by 



deij 



Fv+4 

~f7 



+ 1 



F, 



2 F3 



> 



V + 4 V F, 



v-l rp2 
v+A \ ^v+2 
i?2 



(3.16) 



*from now on we shall drop the argument of the x^-c.d.f. 's, which is always p/A. 
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2 4 6 8 10 12 14 16 18 20 
Fig. 1 - det vs. p/A at w = 3, . . . , 7. 



the lower bound on the r.h.s. following from the inequality Fy^^F^/ F^_^_2 > {v + 2)/(u + 4), first 
proved i n [6]. Accordingly, we conclude that J is non-singular for any finite value of p/X, and there- 
fore eq. (3.14) is unambiguously solved by A*^"-* = J'^G^^^ Note as well that lim^^^_j^g det J" = 0, 



so the invertibility of becomes critical at extremely small values of p/X. By way of example, we 
show in Fig. [l] a plot of det{J') vs. p/X for v = 3, ... ,7. Finally, the inverse of can be readily 
checked to be 



)jk 



F, 

Fv+A 



[{v + 2)5 - 1] 



(^^ + 2: 



Fv+i 

Fv + 4: 



[v6jk - 



1] 



V 



(3.17) 



Let us now come to the analytic structure of the known terms Q 



in) 



owing to the chain rule of differentiation, every single contribution to (except for /^[."'*) 



. We have just explained that 

holds 

a partial derivative di^^ . . . di^Rk{p; At) with some indices ii, . . . , it- On expanding this in terms of 
a and a^, we produce ratios with numerators made of products of derivatives of a and and 
denominators amounting to some power of a. From the rules established by Propositions 2.1 and 

(n) 

2.4, it follows that Gl can be represented in full generality as 



n+1 ki 



Fv+2ki • • • • • Fv+2k„+l 



fcl=0fc2=0 fcn+l=0 

K<n+1 



F 



n+l 



(3.18) 



n+l 



with the coefficient A""""*"^ factored out for later convenience. The subscript prescription K < n + l 
to the nested sum has to be understood as a restricting condition to the possible values taken by 
the sum variables ki, . . . , k^+i. For the sake of conciseness, we shall refer collectively to the ratios 
Fv+2ki • • • • • F^j^2kn+i / ^v~^^ x^-TdXios and the coefficients c^."^^ ^^^^ as c-coefs. We observe that 
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the r.h.s. of eq. (3.18) becomes increasingly populated for large values of n. In the few lowest order 
cases, it expands to 



el" 



A:;00 + '^k;10 rp + ^k;ll 772 



+ 2 , (1) + 

"T ''fc;20 p 



(3.19) 



(2) , J2) , J2) 

fc;000 ^ ''A:;100 ^ H;110 ^2 ^3 



(2) i^t;+4 (2) i^t;+4i^i.+2 (2) ^^^.+6 
''fc;200 p '-k;210 p2 ^fc;300 p 



(3.20) 



Without conditioning the sum to -fC < n + 1, the number of summands would be (^^^i) (cf. eq. (1) 
of [5]). Actually that number is much lower, yet {^^^i) can be taken as a (loose) upper bound to 



it. The intricacy of eq. (3.18) is only fake: upon adding separately the degrees of freedom of all 
the x^-c.d.f. 's at numerator and denominator of each 'x^-iaiio and then subtracting the resulting 



numbers yields 2K. Therefore, eq. (3.18) is just a formal way of representing a linear combination 
of x^^ratios, where each x^~c.d.f. has at least v degrees of freedom and the overall algebraic 
sum of degrees of freedom amounts to 2K = 0,2, .. . ,2(n + 1) (with denominators contributing 
negatively). We stress once more that this analytic structure is a direct consequence of the chain 
rule of differentiation alongside with the results established in Propositions 2.1 and 2.4. 

As for the c-coefs, they do not depend on p and can be only determined by direct calculation. In 
spite of this, their dependence upon the Taylor coefficients of A(e) displays a well defined analytic 
structure. In order to disclose it, we must rely upon the notions of physical and perturbative 
dimensions. 

i) First of all, we assume that A has the physical dimension of length (L), for which we adopt the 
notation [A]l = ■ If we also assume [e]L = then it follows [A]l =iA^L = • • . = [A^^^Jl = ■ 
Similarly, we have = and [Rk]i. = Since = from eq. (ETJ it follows [g^'^\ = L^. 



Hence, eq. (3.18) leads us to conclude 



l''fc;fcl...fc„+iJL 



L". 



(3.21) 



As previously explained, c-coefs can only depend polynomially upon the Taylor coefficients of A(e). 
Eq. (3.21) tells us that such polynomials must be linear combinations of monomials in A and the 
directional components of X^^\ . . . , A^""^-*, each monomial having precisely degree n. 

ii) We define the perturbative dimension of a single monomial as the sum of the perturbative orders 
of its factors. More precisely, we set [A]p = [A^-^^Jp = P^, . . . , [A'-'^^Jp = P". Thus, for instance, 
we have [A^(A-^^)^(Aj"'^^)^]p = P^. We remark that A carries no perturbative dimension, but it does 
increase the physical dimension of the monomials. Since ^^"^ is the result of the expansion at 0(e"), 
it follows 



\ (n) 1 _ pn 

The same result applies clearly to each single contributing monomial. 



(3.22) 
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iii) Several monomials within a given c-coef have the same perturbative structure and numerical 
prefactor and differ only by directional indices, e.g. 3A^(A^^^)^(A2^'')^ and 3A^(A4^-')^(A5^^)^. This 
is a consequence of the index structure of Ajj.,,j^: the products of Kronecker symbols contributing 



to the r.h.s. of eq. (2.6) contract the indices of the Taylor coefficients of A(e) in all possible ways, 



thus generating an increasing number of new aggregate structures at each order of the expansion. 
For instance, monomials within c-coefs belonging to the lowest perturbative orders are grouped 
according to 



0(6^) 



V V 

,(1) A._>^Ma)^2. 



1=1 i=l 



(3.23) 
(3.24) 



i=l 



i=l 



i=l 



(2)^2 



4 = 1 



i=l 



Ci.. = E(^rVAf, Cl:4 = E(^l'^) 



(1)^4. 



(3.25) 



i=l 



i=l 



In view of the above considerations, we conclude that all c-coefs at ©(e") with n > 2 can 
be represented in full generality as linear combinations of all possible products of perturbative 



structures under the constraints imposed by eqs. (3.21) and (3.22), i.e. 

(n) _ (".m) /n(") 

k;ki...kn+i / J Tk\...kn+\ k\ra^ 



(3.26) 



with numerical prefactors 7^"'"^''^^^ and perturbative structures O^'.'^ fulfilling [O^'-'^Jl = and 
[^I'^mlp ~ instance, we have 

oil e {(aI^V, a1^)Ci,C?,Ci:2}, (3.27) 



S {(Afc''')'') (Afc'O^'Cl, A^^^Cl, Cl, A^^^Cl:2, Cl:3, 



ClCl:2, AA[,^''A[i^\ \y^l^C,2, ACl2, AA^^-'Cl, AC1C2} , 



(3.28) 



Before working out the expansion at a given order, one should write down a complete set of 
perturbative structures pertaining to that order, such as eqs. (3.27) and (3.28) illustrate for n = 2, 3. 
A preliminary identification of all suitable structures is indeed particularly useful in order to identify 
groups of terms when high order calculations are performed by means of a computer algebra system 
(CAS), as we shall see in sect. 3.3. 



When calculating c^?^^ ^^^^ , many of the coefficients 7^"'"^^^ 
vanishing ones are subject to 



are found to be zero. The non- 
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Proposition 3.1. For n > 2 the coefficients 7^"'"^^ ^ are constrained by 



-1 fci 



(3.29) 



fci=0fc2=0 fc„+i=0 
K<n+1 



Proof. We first note that if ^ G ^('^p ^)) tlien /i G 2?(r Vp' > p. Therefore, it makes sense 



to consider eq. (3.14) as p — )• oo with fi kept fixed. In particular, we have show ed pr eviously that 
limp_^oo <Jkj = ^kj- Moreover, as /9 — )• oo all the x^~J^atios tend to one, thus eq. (3.14) reduces to 



n+l ki kn n+1 fei fc„ 

'^k = ^fc + X] X/ ■ ■ ■ X] ^fc;fci...fc„+i = ^fc + X/ X/ ■ ■ ■ X]">'fci.'..fc„+i^A;;m 



fci=0A:2=0 fc„+i=0 
K<n+1 



fci=0A:2=0 fc„+i=0 m 
X<n+1 



Z—,' ■ ■ ■ Z-^ 'fci...fc„+i 

fci=0fc2=0 fc„+i=0 
_ft:<n+l 



(3.30) 



(n) (n) 

However, in the same limit Xk — >• Pfc; which entails order by order — )• . Hence, we infer that 



the sum of perturbative structures on the r.h.s. of eq. ( 3.30[) van ishes as p — t- oo. Since in general 



ol^^\fl, fi'^^^ . . . ,fi^^ 7^ 0, we must conclude that eq. (3.29) holds true. 



□ 



Now, the first few orders of the perturbative expansion can be worked out easily with a little 
algebra. Doing the calculations is useful to get acquainted with the general structure examined so 
far. The lowest order has been discussed in sect. 1, so we shall concentrate on the perturbative 



corrections to it. As of now, we shall assume that /i(e) is expanded according to eq. (3.4). 



3.1 Perturbative expansion at O(e^) 
Equations {<?^^'* = 0}^^^^ have the explicit form 

The only term we need to calculate is 



R 



(1) _ dRk 



de 



= Y,^f^djRk{p;K). 

e=0 ^ = 1 



(3.31) 



(3.32) 



Actually, we have calculated the first order partial derivatives of Rk in eq. (3.13). Thus, we have 

(3.33) 



jkj Xj = Sfik 



whence we infer qI^^ = Sfi^- Choosing fi = fJ- yields an important simplification: 



15 



Proposition 3.2. If fi = ji, then C,i 



0. 



Proof. It is sufficient to add side by side all eqs. (3.33) for A; = 1, . . . , u to get 

V 

= 0. 



Ci 

2 



+ V 



F2 



k=l 



Since the quantity in square brackets is strictly positive, we must conclude that C,i = 0. 



(3.34) 



□ 



As can be readily understood, annihilating C,i results in a huge simplification of the higher order 
calculations. Indeed, C,i figures in many perturbative structures contributing to ^^"^ for n > 2. For 



instance, the structure basis of eq. (3.27) is reduced to only two elements in place of four, while the 



one of eq. (3.28) is reduced to six elements in place of twelve. 



3.2 Perturbative expansion at O(e^) 

The subleading correction A*^^) is obtained from the equations 



(3.35) 



Most of the contributions to the three terms on the r.h.s. are calculated smoothly at this point. 
For instance, we have 



1 

I L 



(1)a2 



A F^ 



A F^ 



(3.36) 



A^^^i? + \R 



(2) 



2 de2 



e=0 



E^^^-^f + ^ E A5.;u«%,,i?.(p;A,). 



(3.37) 



iii2=i 



To keep things general, we make no assumptions on fi here. Accordingly, we keep all x^-vsXios 
proportional to powers of Ci- The evaluation of the second derivatives dj-^j^RkiP', At) requires a few 
pages of tedious algebraic work, which we cannot detail for obvious reasons. The upshot is given 
by 



i=i 



(2) 



1 

"SA 



Ci' + 2Ci;2 + 4A1'^Ci+8(aW)2 



F, 



+6 



F. 



+ 



1 



8A 



3C? + 2Ci:2 + 4A«Ci 



Fv+aFv+2 



(-2 p3 



+2 



F^ 4A F^ 



+ ^[C.2 + 2(A, )]— 



(3.38) 



(2) 

We notice that the four structures listed in eq. (3.27) all contribute to the c-coefs c^.^^ 



lfc2fc3' 
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77.3 




172 






F2 

V 


^ V 




V 




-1 








1 





^ (1) * 


1 

^2 


1 
2 











a 


1 

^8 


CO 1 00 


1 

~4 








Cl:2 


1 

^4 


1 

4 





1 

2 


1 

~2 



Table 1 - Coefficients Ik'ti^k^- 



In Table 



f2 rn) 

we collect the coefficients 'Jf^-^'j^^j^^- Instead of naming rows and columns respectively 
according to tne values of m and the triples {ki, k2, ks), for the sake of readability we identify each 
table entry by the perturbative-structure and the x^^ratio it refers to. This way of tabulating the 
coefficients becomes particularly informative at higher orders. We ffiially observe that adding the 
entries of each table row yields zero, in accordance with eq. (3.29). 



3.3 Perturbative expansion at higher orders 

Paper-and-pencil calculations become prohibitively expensive at higher orders. Fortunately, it is 
not difficult to work out the algebra with the assistance of a CAS. For the reader's convenience, 
in the appendix we attach prototype Maple"'''^ procedures, which facilitate the task. The code is 
split into three blocks, which we shortly review. 

The first code block (A.l) contains a procedure DeltaO , which computes the coefficient Aj^,..j„. 
The procedure argument is assumed to be a list of nonnegint items; alternatively the procedure 
returns unevaluated. The input list is ffist sorted in ascending order, then the multiplicity set is 



identified. The procedure computes the r.h.s. of eq. (2.2) and returns its numerical value 



The second code block (A. 2) performs the algebraic work related to the perturbative expansion 



of eq. (1.1). Before submitting it to evaluation, the user is assumed to assign a nonnegint variable v 
representing the number of dimensions, and a nonnegint variable n < 4 representing the highest 
perturbative order processed by the program. The code block starts with a pair of procedures. 



DerAlphaO and DerAlphakO, which encode respectively eqs. ( |2.24 ) and (2.25). Afterwards, it 



performs a Taylor expansion of the r.h.s. of eq. (1.1) up to O(e^). Taylor coefficients are stored 



within an indexable object h[j ,k] , the indices j and k representing respectively the perturbative 
order and the physical direction. At this stage, h[j,k] holds a sum of potentially many terms. 
The summands contain derivatives of a and ak, which are purely symbolic objects. Their evalua- 
tion requires sequences of prescriptions, stored within the the variables COA, COAk,. . . , C4A, C4Ak. 
Algebraic simplifications are performed in the last few lines, where partial results are stored within 
indexable objects hOO,. . . , hlO, so as to allow for an offline oversight of the single steps. 



The third code block (A. 3) illustrates in a specific case a numerical technique which we have 

) 

n + l ■ 



devised for the determination of the coefficients 7^."'"^^ . The code processes the coefficients 
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Pv+8 


-Pi)+6^D+2 


P2 


Fv+6 


Fy^4Fy^2 




T7'2 

-^^11+2 








p2 

^ V 




p2 

^ V 


F. 






-1 








2 





-1 







1 

^4 





1 

4 


1 
2 


1 

~2 








Cl:3 


1 

^6 


1 

6 





1 
2 


1 

~2 


1 


1 

2 


AAi > 











-2 





2 





aaWc2 











1 

~2 


1 
2 








AC12 











1 

~2 


1 

2 


1 


-1 



Table 2 - Coefficients l^k^kl^k^ki- assume in tliis case Ci — 0. 



corresponding to n = 3 and {ki, k2, k^, k^) = (3,0,0,0), i.e. those entering the c-coef multiplying 
the x^^ratio F^^q/F^. It also assumes fi = p,, which reduces the basis of perturbative structures to 



.(1)^3 x{l) 



(l)x(2) 



(3.39) 



The algebraic sum pointed to by hlO [3,k] at the end of the second code block has no knowledge 
of these structures. In order to identify them within hl0[3,k] , we need to group terms properly. 
Instead of proceeding at an algebraic level, which turns out to be computationally demanding, we 



adopt a numerical approach, based on the use of eq. (3.26) as a square linear system fulfilled by the 



coefficients 7l'^'"l^ . Having subtracted from hl0[3,k] all contributions figuring in eq. (|3.12|), we 
extract from it all terms proportional to F^^q/F^, whose sum amounts to — c 



(3) 

fc;3000' 



Then, for each 



0), from which we compute O 



(3) 

fc;l' ■ • • ' 



k we assign A^^^ and A*-^^ random values (chosen so that 

(3) (3) (3) 

O^.g and c^.3ooo- "^^^ random matrix O^.^ thus obtained is non-singular, so eq. (|3.12|) can be solved 



with respect to 73000^- The solution must be independent of the random numbers extracted. This 
represents a strong hint of goodness of our determination, but a real check consists of an algebraic 

(3) 

comparison between the reconstructed coefficient c^.jqqq and the one extracted from hl0[3,k] . In 



Tables |2|and|3| we report the coefficients Tfcj^fc2fc3fc4 '^^^ ^fci'fc2fe3fc4fc5 ^'^'^^-'^ the assumption 



0. 



4 Properties of the first few perturbative coefficients 

So far we have focused on formal aspects of the perturbative expansion with the aim of proving its 
theoretical and computational feasibility. When it comes to establish the level of accuracy reached 
by approximating the reconstruction with a finite number of contributions, we are soon led to 
investigate the analytic properties of the perturbative coefficients of A. 
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Hereinafter, we shall consider perturbative series truncated at different orders, for which we find it 
worthwhile introducing the notation 



a1"U^A«. (4.1) 

1=0 



Fig. [2] shows an illustrative example of perturbative reconstructions at f = 4, which we shall focus 
on in this and next section. Plots have been produced as follows. First of all, in order to probe 
perturbation theory on an eigenvalue spectrum characterized by a relatively large condition number, 
we have chosen Aex = {0.1,0.3,0.8,2.2} as the full eigenvalue spectrum to reconstruct. By means 
of numerical techniques detailed in [1], we have computed fi = Tp ■ Aex for several values of p. In 
correspondence with each pair {p, p) we have finally reconstructed the eigenvalue spectrum up to 
the fourth perturbative order, having chosen in all cases p = p. We notice from the plots that 
the error made upon truncating the expansion at a given order increases in general at lower values 
of p. Moreover, the error is larger for eigenvalues at the edges of the spectrum and milder in the 
center of it. Remarkably, the convergence pattern of the lower half of the spectrum is radically 
different from the upper half. In the former case the perturbative series seems indeed to converge 
with alternate signs, whereas in the latter it displays an almost monotonic character. 

In order to explain the observed behavior, let us first concentrate on the leading contribution A. 



If {p*,p) fulfills the constraint < p < p*{v + 2)~^, a solution A to eq. (1.3) exists for all pairs 
(p, p) with p > p* . In this sense we can consider A as a function of p at fixed p. The analytic form 
of Tp tells us that limp^oo X = p- Now, we know for sure (modulo conjectures on the structure of 
&b) that p < p*{v + 2)-^ provided p £V {TpJ'). Since pi < p < Pv, continuity we conclude 
that 

3p: Tp-\pi)<\<Tp\p,) yp>p. (4.2) 
At sufficiently large p the above inequality holds true of cou rse w ith pi and py respectively replaced 



by Pi and pi+i, where i is such that pi < p < pi+i- Eq. (4.2) does not tell where A is placed in 
relation to the full eigenvalue spectrum. To find such an estimate, we can resort to eq. (2.12) of [1]. 
Based on arguments which are completely analogous to those used in there, we arrive easily at 

p M{v,v + 3/2,p/{2Xi)) ^ ^^^^^ ^ p M(l,5/2,/>/(2A.)) ^^^^ 



2v + lM{v,v + l/2,p/{2\i)) - ' - 3M(l,3/2,/9/(2A„)) 
with M(a, 6, z) denoting a Kummer function, viz. 



ri=0 



As a consequence of the asymptotic limit 7^ (A) ^——^ A and (see e.g. chap. 13 of [7]) 

p M{v,v + 3/2,p/{2X,)) p^oo . 



2v + lM{v,v + l/2,p/{2Xi)) 
pM{l,5/2,p/{2Xy)) p^ 



3M(l,3/2,/5/(2A,)) 



A. , (4.6) 
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,(1) 



10 12 14 16 18 20 22 24 



,(2) 



.(4) 



> P/A 



0.36 




10 12 14 16 18 20 22 24 



.(2) 



.(3) 



,(4) 



p/A 




0.65 



,(1) 



10 12 14 16 18 20 22 24 



,(2) 



.(4) 



> P/A 




1.60 



> P/A 



Fig. 2 - An example of perturbative eigenvalue reconstruction at w = 4. The unconstrained 
spectrum Aox = {0.1,0.3,0.8,2.2} is represented by black dashed lines. The solid lines corre- 
spond to the four levels of approximation obtained by truncating the perturbative series at the 
first to fourth order. 



we conclude that if p is sufficiently large, then Ai < A < A^,. As intuitively expected, non-linear 
effects are mitigated in the region of weak truncation. The situation is qualitatively depicted 
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Fig. 3 - A schematic diagram showing the action of the operators ^ (black dashed hnes) 
and 7p~^ (green thick dashed hne) when p. = fl and ^ G 'D{Tp^). Ifl<i<t)~lis such that 
^ P- l£ Mi+ii then at sufficiently large p we also have Ai < A < A^+i. 



m 



Fig. [31 



Let us now consider the first few corrections to A. As far as we are concerned with their 



numerical computation, we can hmit ourselves to invert eqs. (3.18) one by one by means of a linear 



solver. On the other hand, the matrix structure of J^~^ (whose off-diagonal entries are all the 
same) allows us to perform the inversion in closed form. If the expression obtained as a result has 
too many contributions, we shall hardly find it better. Certainly, this is not the case with the first 
few perturbative corrections, of which we want to estimate the range of variation. 



Since O 



is either index-free or dependent upon k via monomials (A^ 



with 



iimi + . . . + irmr < n, the only algebraic ingredients we need for the analytic inversion of eqs. (3.18) 
are the sums 



jk 



2V- 



V = {v + 2)- 



v+i 



k=l 



^ V 



(4.7) 



(.ir)\mr 



Cii:mi...ir:mr'^ 



-1 f F]j+4 



' v+2 



where Cn:mi...v:m, = EJ=i(Af ^)™^ • • • (Af ^)'^'- is defined in perfect analogy with eqs. (|3.23|)-( 



(4.. 



We are thus ready to work out the algebra. In first place, a straightforward calculation yiek 



3.25). 



A 



(1) 



whence we infer that 
\i = Ji = 



sign (A^^^) = sign {fij, - /i) 



(4.9) 
(4.10) 
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We thus recognize that the first perturbative correction to the leading term A is positive for the 
upper part of the spectrum and negative for the lower one. 

The algebraic evaluation of the second perturbative correction to A is as easy as the first one, 
yet estimating its range of variation is somewhat more difficult. Let us set in this case jl = fi 
from the beginning. Under this assumption, we see that the r.h.s. of eq. (3.38) reduces to a linear 
combination of four terms, namely 



m=l 
1 

~4A 



(2,m) 
k 



Cl:2+4(A«)^ 



v+6 



^[Cl:2 + 2(A, )] — 



On inverting eq. (3.38) we find 



.(2) 



m=l 



g(2,4) 



Cl:2 



4A 



F2 



Cl:2 Fv+2 

' 2~X 



(4.11) 
(4.12) 
(4.13) 

(4.14) 



>(2,1) 



,(2,2) 



,(2,3) 



A(2,4) 
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Hence, we have 
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Since F„+e < F^,+4, the first contribution to the r.h.s. is certainly positive. As for the second one, 
we define 

^ Fv+6 
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Plots reported in Fig. [4] suggest H > and consequently A^^^ > 0. Unfortunately, we lack at the 
moment an analytic proof of such inequality. Nevertheless, if the numerical evidence is correct, we 
can conclude that all the eigenvalues receive a positive contribution from the second perturbative 
correction, which together with eq. (4.10) explains qualitatively why the lower part of the spectrum 
seems to converge with alternate signs, whereas the upper part has an almost monotonic character. 
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4 8 12 16 20 24 28 32 

Fig. 4 - Positivity of the linear combination S of x^^ratios. 



5 Perturbative estimators vs. the iterative one 



In the last part of the paper we introduce a source of statistical uncertainty. So far we have studied 
the eigenvalue reconstruction under the hypothesis that fi = Tp ■ X represents the exact truncated 
counterpart of some A G M''. This is rather unusual in most applications, where // is not the result 
of an exact truncation, but is instead estimated from a normal population Vn = {^^'•'^^IfcLi of 
finite size N, distributed with S = diag(A). The sample estimate we consider here is performed 
as follows: a certain subset of M < elements of Vn falls within Bv{p), with the fraction M/N 
fulfilling limAr_s.oo M/N = a{p; A). From this subset we measure ©g via the classical estimator 



N 



k=l 



Xi 



^j:xf).i[xWG^.(p)]. 



(5.1) 
(5.2) 



k=l 



The eigenvalue spectrum /i of represents our definition of the sample estimate of /u, which we 
use as an input parameter for the reconstruction of A. Of course, fi is interpreted as the realization 
of a stochastic variable in sample space. It thus makes sense to pose the question of what the 
statistical properties of the stochastic variable A = • fi and its perturbative approximations 
AW (A; = 1,...,4) are. 

Finding analytic relations between the expectation value in sample space of polynomial functions 
of fi and analogous functions of A is a difficult task, since we dispose of no analytic representation 
of the reconstruction operator r~^. This goes beyond the aims of the present paper, so we limit 
ourselves to perform a simulative study in the specific case where Vn distributes according to 
S = diag(Aex), with Aex introduced in the previous section. In our study, we have chosen = 
200, 250, . . . , 2000; for each value of A^, we have generated about 5000 normal populations; for 
each of them, we have then considered Euclidean balls with p = 4.0, 6.0, . . . , 40.0 and for each 
pair {p, N) we have finally measured bias and variance of A and A(^). As an example, in Fig. pi we 
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report our results at p = 6.0 (corresponding to a weak truncation with a(6.0; Aex) — 0.844|^ From 
the plots on the left we notice that 

i) the bias of A^^^ is weakly sensitive to N for all fc's; it converges asymptotically to the intrinsic 
perturbative bias corresponding to /U = limjv_s.oo with finite size corrections proportional 
to 1/iV; 

ii) A is slightly biased at finite A'" and asymptotically unbiased; convergence is again reached 
linearly in 1/N. 



Similarly, from the plots on the right we observe that 



iii) all variances vanish linearly in 1/N; 

iv) the variance of the higher eigenvalues increases at fixed A'^ as we add perturbative corrections; 

v) by contrast, the variance of the lower eigenvalues decreases at fixed as we add perturbative 
corrections; 

vi) the iterative estimator has a higher variance than all its perturbative approximations at the 
top of the spectrum and a lower one at the bottom of it. 



The variance plots illustrate the potential usefulness of the perturbative estimators. The recon- 
struction of the upper part of the spectrum achieved from the iterative algorithm is rather noisy 
when N is moderately small (e.g. in this case N < 400). The perturbative estimators allow to 
control the variance, the price to pay being the introduction of an asymptotic non-vanishing bias. 
Depending on the specific context, there will be an optimal choice for the order of the perturbative 
approximation, which guarantees acceptable values of both bias and variance. 

Results are qualitatively similar for the other simulated values of p: the asymptotic biases and 
the slopes of the variances decrease as p increases, as naively expected. 

A noticeable outcome of our simulations is inferred upon relating the variances of the recon- 
structed eigenvalues to those of the truncated ones. To this aim, we observe that since A = ■ p 
is a vector equation, each component of the reconstructed spectrum Aj = Xi{pi, ...,//«) depends 
upon all the components of p. Hence, it follows that var(Aj) is itself a function of all the compo- 



nents of p. Nevertheless, if eq. (3.8) holds true, Aj depends weakly on pj^ for k ^ i. Therefore, it 
makes sense to look at how var(Aj) relates to var(/ij). An example of such dependence is shown 
in Fig. [6]for i = 1,4, corresponding respectively to the lowest and highest components of Aex- We 



first note that points displace along definite curves, in accordance with the conjecture of eq. (3.8). 
Then we note that the variances are linearly related, except for weak nonlinear effects observed at 
var(/i4) ^ 1.0 X 10-2. A gain, we observe that the variance of the iterative estimator of the lowest 
eigenvalue is minimal and that of the highest one is maximal. What is most remarkable is the slopes 
observed for the highest eigenvalue A4. On comparing the scales of the x- and y-axis, we recognize 
that a huge inflation of the variance occurs as a result of applying to p. Numerical simulations 

^statistical errors of the sample estimate of the variances have been computed according to the general formula 
for the standard error se(var(AA;)) = var(Afe) ■ [2/{N — 1) -\-k/N], where var(Afe) is the sample estimate of var(Afe) and 
k is the sample excess kurtosis of the distribution of var(Afc). 
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signal the existence of such amphfication phenomena, yet they are not able to unveil the underlying 
mechanism. One should be anyway aware that in practical situations the exact reconstruction of 
the highest eigenvalue may be critical. In such cases, the adoption of perturbative estimators in 
place of the iterative one may represent a good way out. 

6 Conclusions 

In this paper we have explored a perturbative approach to the reconstruction of a normal covariance 
matrix S from a spherically truncated counterpart ©g. We recall that S and &b commute, so the 
reconstruction problem concerns only the eigenvalue spectra. Having preliminarily collected all the 
ingredients needed for the implementation of the perturbative expansion, we have detailed both its 
general features and the practical aspects related to the calculation of the perturbative coefficients. 
In the paper we provide formulae for the reconstruction of the eigenvalues of S up to the fourth 
perturbative order as well as a prototype Maple"^'^ code to further improve the approximation. 

Prom a theoretical point of view, the perturbative method is meant to complement the fixed- 
point iterative algorithm proposed by us in ref. [Ij in cases where the latter becomes inefficient. 
Such cases occur when p is comparable or less than the lowest eigenvalue of S and/or the number 
of dimensions v is very large. In both cases, the inefficiency consists in a slow convergence speed. 
A second limit of the iterative algorithm shows up when the eigenvalue reconstruction is performed 
from statistically poor sample estimates of the truncated covariance spectrum. Specifically, we 
have shown that a statistical uncertainty on the truncated eigenvalues is inflated by the application 
of the reconstruction operator, thus producing large fluctuations of the higher components of the 
reconstructed spectrum. Perturbation theory offers the possibility to control the variance and 
stabilize the reconstruction by properly choosing the approximation order. The price to pay upon 
replacing the iterative estimator with its perturbative approximations is the introduction of an 
asymptotic non-vanishing bias. Of course, it is possible to adopt mixed approaches, where the 
lower part of the spectrum is reconstructed via the iterative estimator while the upper part is 
obtained from a perturbative one. 

We conclude by recalling that both approaches developed in this paper and ref. |lj are based on 
some conjectured correlation inequalities over Euclidean balls, which have been extensively tested 
via numerical simulations and will be further analyzed elsewhere. 
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as 
> 
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Fig. 6 - Variances of the iterative and perturbative estimators of the lowest and highest 
reconstructed eigenvalue vs. the variance of the corresponding truncated eigenvalue. Normal 
populations have been generated with S — diag(Aox)- 
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Appendix A Maple^m code 
A.l Code block 1: the coefficient Aj^ 



# Delta coefficient 
# 

Delta := procO 

local Sort Args , V , ActCtr , NxtCtr , Res , k : 
for k from 1 to _npassed do 

if not type (_passed[k] , 'nonnegint ' ) then 

return 'procname (_passed) ' : 
end if : 
end do : 

SortArgs := sort ( [_passed [1 . . _npassed] ] ) : 
V := Vector(_npassed) : 
V[l] := 1: 
ActCtr := 1: 
NxtCtr := 2: 

for k from 1 to (_npassed-l) do 

if SortArgs [NxtCtr] = SortArgs [ActCtr] then 

V [ActCtr] := V [ActCtr] +1: 
else 

ActCtr := NxtCtr: 

V [ActCtr] := 1: 
end if : 

NxtCtr := NxtCtr+1: 
end do : 
Res := 1: 

for k from 1 to _npassed do 

Res := Res*(doublefactorial(2*V[k]-l)) : 
end do : 
return Res : 
end: 

A. 2 Code block 2: perturbative expansion of eq. (|l 

# Nested sequence 
# 

NestSeq := proc (TheEq.v : :nonnegint , niter : :nonnegint) 
if niter = then 

eval(TheEq) : 
else 

seqCeval (NestSeq (TheEq.v, niter- 1) ) , 
catCr', niter) = l..v): 

end if 
end proc : 

# Derivatives of Gaussian Integrals 
# 

DerAlpha := procO 
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global V, Delta: 
local m.Factl ,Fact2 : 
m := _npassed: 

Factl := Delta(_passed[l. ._iipassed])/(2*l[0])~m: 
Fact2 := add((-l)~(m-j)*binomial(m, j)*F[v+2*j] , j=0. .m) : 
return Factl*Fact2: 
end proc: 

DerAlphai := procO 
global V, Delta: 
local m.Factl ,Fact2: 
m := _npassed-l: 

Factl := Delta(_passed[l. ._npassed])/(2*l[0])~m: 
Fact2 := add((-l)~(m-j)*binomial(m, j)*F[v+2*(j+l)] , j=0. .m) : 
return Factl*Fact2: 
end proc: 

# Function arguments 
# 

lam : = Vector (v) : 
for k from 1 to v do 

lam[k] := add(l [j ,k] *epsilon~ j , j=0. .n) : 
end do: 

lam := seqClam [k] ,k=l . . v) : 
lamO := seq(l [0,k] ,k=l . .v) : 

# Integral ratio 
# 

R := proc(j) 
global lam: 

return alpha [j] (lam) /alpha(lam) : 
end: 

# Taylor expansion of the map 
# 

for j from 1 to n do 
for k from 1 to v do 

Rk := convert(taylor(R(k) ,epsilon=0,n+l) ,polynom) : 
h[j,k] := expandCcoeff (lam[k]*Rk,epsilon, j)) : 
end do : 
end do: 

# Evaluation conditions 
# 

COA := alpha(lamO)=F[v] : 

COAk := seq(alpha[k] (lamO)=F[v+2] ,k=l. .v) : 

CIA := NestSeq(D[rl] (alpha) (lamO)=DerAlpha(rl) ,v,l) : 

ClAk := NestSeq(D[rl] (alpha[r2] ) (lamO)=DerAlphak(rl,r2) ,v,2) : 

C2A := NestSeq(D[rl,r2] (alpha) (lamO)=DerAlpha(rl ,r2) ,v,2) : 

C2Ak := NestSeq(D[rl,r2] (alpha [r3] ) (lamO)=DerAlphak(rl,r2,r3) ,v,3) : 
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C3A := NestSeq(D[rl,r2,r3] (alpha) (lamO)=DerAlpha(rl,r2,r3) ,v,3) : 

C3Ak := NestSeq(D[rl,r2,r3] (alpha[r4]) (lamO)=DerAlphak(rl,r2,r3,r4) ,v,4) : 



C4A := NestSeq(D[rl,r2,r3,r4] (alpha) (lamO)=DerAlpha(rl ,r2,r3,r4) ,v,4) : 

C4Ak := NestSeq(D[rl,r2,r3,r4] (alpha [r5] ) (lamO)=DerAlphak(rl,r2,r3,r4,r5) ,v,5) : 

CArgO := seq(l [0,k] =1 [0] ,k=l . . v) : 

# Evaluations 
# 



for j from 1 
for k from 
hOO [j ,k] 
hOl [j ,k] 
h02[j,k] 
h03[j,k] 
h04 [j ,k] 
h05[j ,k] 
h06[j ,k] 
h07[j,k] 
h08[j,k] 
h09[j,k] 
hlO[j,k] 
end do : 

end do: 



to n do 

1 to V do 

expand (eval(h[j 
expand (evaKhOO 
expEind ( e val (hO 1 
expand (eval (h02 
expand (eval(h03 
expand (eval (h04 
expand (eval (h05 
expand ( eval (h06 
expand (eval (h07 
expand (eval (h08 
expand ( e val (h09 



,k] , [CIA])) : 
[j,k] , [C2A])) : 
[j,k],[C3A])): 
[j.k] ,[C4A])): 
[j,k] , [ClAk])) 
[C2Ak])) 
[C3Ak])) 
[C4Ak])) 
[COA])): 
[COAk])) : 
[CArgO] ) ) 



[j,k]. 
[j,k]. 
[j .k] : 

[j,k]. 
[j,k] , 



A. 3 Code block 3: extraction of 7^"'^^ 



with (Linear Algebra) : 
with(RandomTools) : 

V := 6: 

# 0-structure matrix 
# 



zetal := add(l[l,k] ,k=l. .v) : 
zeta2 := add(l [2 ,k] ,k=l . . v) : 
zetall := add(l [1 ,k] -2,k=l . . v) : 
zetal2 := add(l[l,k]*l[2,k] ,k=l. 
zetalll := add(l[l,k]-3,k=l. .v) : 



.v): 



SSmatrix := Matrix(v,v) : 



for k from 1 to v do 



SSmatrix [k, 1] 
S3matrix[k,2] 
S3matrix[k,3] 
S3matrix [k,4] 
SSmatrix [k,5] 
SSmatrix [k, 6] 
end do : 



= l[l,k]-3: 

= expandd [1 ,k] *zetall) : 

= expand (zetalll) : 

= l[0]*l[l,k]*l[2,k] : 

= expandd [0] *1 [1, k] *zeta2) : 

= expandd [0] *zetal2) : 
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# Jacobian matrix 
# 

Jmatrix := Matrix(v,v) : 

for kl from 1 to v do 

for k2 from 1 to v do 

Jmatrix [kl,k2] := (l/2)*(Delta(kl,k2)*F[v+4]/F[v] - F[v+2] -2/F[v] "2) : 

end do: 
end do: 

# Terms to be removed by hand 
# 

V := Vector (v) : 

for j from 1 to v do 

V[j] := 0: 

for k from 1 to v do: 

V[j] := V[j] + Jmatrix [j,k]*l [3, k] : 
end do: 
end do: 

# Randomized Linear system 
# 

C := Vector (v) : 

for j from 1 to v do 

lincondO : = 1 [0] =Generate (float (range=0 . . 1 , 'method=\mif orm' ) ) : 
linvalsl := seqCGenerate (float (range=0 . . 1 , 'method=unif orm' ) ) ,m=l . . v-1) : 
lincondl := seq(l [1 ,m] =linvalsl [m] ,m=l . . v-1) : 
lincondl := lincondl , 1 [1 ,v] =-add(linvalsl [m] ,m=l .. v-1) : 

lincond2 := seq(l [2, m]=Generate (float (range = . . 1 , 'method=im.if orm' ) ) ,m=l . . v) : 
for m from 1 to v do 

SSmatrix [j ,m] : = eval (SSmatrix [j ,m] , [lincondO , lincondl , lincond2] ) : 
end do: 

r := expand (hlO [3. j] - V[j]): 
s := expand(eval(coeff (r,F[v+6] ) ,F[v+2]=0)) : 
C[j] := eval ((1[0] "2) *F[v]*s, [lincondO, lincondl, lincond2] ) : 
end do: 

# (-1) x Gamma coefficients 
# 

Gcoefs := LinearSolve(S3matrix,C) : 
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